Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db3) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")
6.3.1.1.1 Number of samples used
[1] 27
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.2 Richness

betadisper(beta_div_richness_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.971
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.973
Hot_dry   0.97302        
adonis2(formula=beta_div_richness_wild$S ~ Population, data=wild[labels(beta_div_richness_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.542719 0.2095041 6.625717 0.001
Residual 25 5.820951 0.7904959 NA NA
Total 26 7.363669 1.0000000 NA NA

6.3.1.3 Neutral

betadisper(beta_div_neutral_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.938
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.937
Hot_dry   0.94763        
adonis2(formula=beta_div_neutral_wild$S ~ Population, data=wild[labels(beta_div_neutral_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.918266 0.2608511 8.822682 0.001
Residual 25 5.435610 0.7391489 NA NA
Total 26 7.353876 1.0000000 NA NA

6.3.1.4 Phylogenetic

betadisper(beta_div_phylo_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999   0.11
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.11
Hot_dry   0.12705        
adonis2(formula=beta_div_phylo_wild$S ~ Population, data=wild[labels(beta_div_phylo_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.3218613 0.2162815 6.899207 0.001
Residual 25 1.1662981 0.7837185 NA NA
Total 26 1.4881594 1.0000000 NA NA

6.3.1.5 Functional

betadisper(beta_div_func_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq     F N.Perm Pr(>F)
Groups     1 0.019387 0.019387 1.653    999  0.222
Residuals 25 0.293200 0.011728                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.213
Hot_dry   0.21033        
adonis2(formula=beta_div_func_wild$S ~ Population, data=wild[labels(beta_div_func_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.0831048 0.1680538 5.05002 0.048
Residual 25 0.4114083 0.8319462 NA NA
Total 26 0.4945131 1.0000000 NA NA
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")
6.3.2.0.1 Number of samples used
[1] 27
beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.1 Richness

betadisper(beta_div_richness_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.11796 0.117959 12.963    999  0.001 ***
Residuals 25 0.22748 0.009099                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.004
Hot_dry  0.0013711        
adonis2(formula=beta_div_richness_accli$S ~ Population, data=accli[labels(beta_div_richness_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.639807 0.179834 5.481634 0.001
Residual 25 7.478640 0.820166 NA NA
Total 26 9.118447 1.000000 NA NA

6.3.2.2 Neutral

betadisper(beta_div_neutral_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07844 0.078443 5.2384    999  0.028 *
Residuals 25 0.37437 0.014975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.025
Hot_dry  0.030815        
adonis2(formula=beta_div_neutral_accli$S ~ Population, data=accli[labels(beta_div_neutral_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.947003 0.2306127 7.493387 0.001
Residual 25 6.495736 0.7693873 NA NA
Total 26 8.442739 1.0000000 NA NA

6.3.2.3 Phylogenetic

betadisper(beta_div_phylo_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.06739 0.067395 2.9532    999  0.089 .
Residuals 25 0.57052 0.022821                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.083
Hot_dry  0.098068        
adonis2(formula=beta_div_phylo_accli$S ~ Population, data=accli[labels(beta_div_phylo_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.2441653 0.1224638 3.488854 0.021
Residual 25 1.7496100 0.8775362 NA NA
Total 26 1.9937754 1.0000000 NA NA

6.3.2.4 Functional

betadisper(beta_div_func_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)
Groups     1 0.02351 0.023513 0.635    999  0.442
Residuals 25 0.92569 0.037028                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.429
Hot_dry   0.43303        
adonis2(formula=beta_div_func_accli$S ~ Population, data=accli[labels(beta_div_func_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.0279416 0.024809 0.6360037 0.452
Residual 25 1.0983283 0.975191 NA NA
Total 26 1.1262699 1.000000 NA NA
beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")
6.3.3.0.1 Number of samples used
[1] 54
beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)
6.3.3.0.2 Richness
betadisper(beta_div_richness_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05014 0.050145 6.2252    999  0.022 *
Residuals 52 0.41886 0.008055                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                          0.02
1_Acclimation 0.015808              
adonis2(formula=beta_div_richness_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_neutral_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.799791 0.222218 4.761789 0.001
Residual 50 13.299591 0.777782 NA NA
Total 53 17.099381 1.000000 NA NA
6.3.3.0.3 Neutral
betadisper(beta_div_neutral_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0199 0.0199035 2.1213    999  0.172
Residuals 52 0.4879 0.0093827                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.172
1_Acclimation 0.15128              
adonis2(formula=beta_div_neutral_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_neutral_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 4.770321 0.2856195 6.663569 0.001
Residual 50 11.931346 0.7143805 NA NA
Total 53 16.701667 1.0000000 NA NA
6.3.3.0.4 Phylogenetic
betadisper(beta_div_phylo_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01334 0.013340 0.6524    999  0.441
Residuals 52 1.06332 0.020449                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.451
1_Acclimation 0.42294              
adonis2(formula=beta_div_phylo_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_phylo_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.855070 0.2267502 4.887385 0.001
Residual 50 2.915908 0.7732498 NA NA
Total 53 3.770978 1.0000000 NA NA
6.3.3.0.5 Functional
betadisper(beta_div_func_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01264 0.012640 0.4951    999  0.511
Residuals 52 1.32764 0.025532                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              0_Wild 1_Acclimation
0_Wild                       0.521
1_Acclimation 0.4848              
adonis2(formula=beta_div_func_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_func_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.1558147 0.0935514 1.720109 0.324
Residual 50 1.5097366 0.9064486 NA NA
Total 53 1.6655513 1.0000000 NA NA
beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Do the antibiotics work?

6.3.4.1 Antibiotics

treat1 <- meta  %>%
  filter(time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat1))]
identical(sort(colnames(treat1.counts)),sort(as.character(rownames(treat1))))

treat1_nmds <- sample_metadata %>%
  filter(time_point == "2_Antibiotics")
6.3.4.1.1 Number of samples used
[1] 23
beta_div_richness_treat1<-hillpair(data=treat1.counts, q=0)
beta_div_neutral_treat1<-hillpair(data=treat1.counts, q=1)
beta_div_phylo_treat1<-hillpair(data=treat1.counts, q=1, tree=genome_tree)
beta_div_func_treat1<-hillpair(data=treat1.counts, q=1, dist=dist)
6.3.4.1.2 Richness
betadisper(beta_div_richness_treat1$S, treat1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.015319 0.0153186 6.8764    999  0.014 *
Residuals 21 0.046782 0.0022277                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.015
Hot_dry  0.015919        
adonis2(formula=beta_div_richness_treat1$S ~ Population, data=treat1[labels(beta_div_neutral_treat1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.356644 0.1527052 3.784762 0.001
Residual 21 7.527429 0.8472948 NA NA
Total 22 8.884073 1.0000000 NA NA
6.3.4.1.3 Neutral
betadisper(beta_div_neutral_treat1$S, treat1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.030536 0.0305358 3.8593    999  0.054 .
Residuals 21 0.166158 0.0079123                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.068
Hot_dry  0.062842        
adonis2(formula=beta_div_neutral_treat1$S ~ Population, data=treat1[labels(beta_div_neutral_treat1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.785669 0.2085055 5.532084 0.001
Residual 21 6.778468 0.7914945 NA NA
Total 22 8.564137 1.0000000 NA NA
6.3.4.1.4 Phylogenetic
betadisper(beta_div_phylo_treat1$S, treat1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.012041 0.012041 0.9898    999  0.372
Residuals 21 0.255459 0.012165                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.366
Hot_dry   0.33111        
adonis2(formula=beta_div_phylo_treat1$S ~ Population, data=treat1[labels(beta_div_phylo_treat1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.8963254 0.1888758 4.889993 0.001
Residual 21 3.8492558 0.8111242 NA NA
Total 22 4.7455811 1.0000000 NA NA
6.3.4.1.5 Functional
betadisper(beta_div_func_treat1$S, treat1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01969 0.019691 0.4738    999  0.528
Residuals 21 0.87274 0.041559                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.529
Hot_dry   0.49877        
adonis2(formula=beta_div_func_treat1$S ~ Population, data=treat1[labels(beta_div_func_treat1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.0246208 0.0133549 0.2842492 0.64
Residual 21 1.8189576 0.9866451 NA NA
Total 22 1.8435784 1.0000000 NA NA
beta_richness_nmds_treat1 <- beta_div_richness_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat1 <- beta_div_neutral_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat1 <- beta_div_phylo_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat1 <- beta_div_func_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = join_by(sample == Tube_code))

6.3.4.2 Acclimation vs antibiotics

treat <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
6.3.4.2.1 Number of samples used
[1] 50
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)
6.3.4.2.2 Richness
betadisper(beta_div_richness_treat$S, treat$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.025318 0.0253178 6.021    999  0.011 *
Residuals 48 0.201837 0.0042049                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.015
2_Antibiotics      0.017817              
adonis2(formula=beta_div_richness_treat$S ~ time_point*Population, data=treat[labels(beta_div_neutral_treat$S),], permutations=999,strata=treat$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 4.885035 0.2455889 4.991572 0.001
Residual 46 15.006068 0.7544111 NA NA
Total 49 19.891103 1.0000000 NA NA
6.3.4.2.3 Neutral
betadisper(beta_div_neutral_treat$S, treat$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.039587 0.039587 6.8387    999  0.008 **
Residuals 48 0.277854 0.005789                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.008
2_Antibiotics      0.011886              
adonis2(formula=beta_div_neutral_treat$S ~ time_point*Population, data=treat[labels(beta_div_neutral_treat$S),], permutations=999,strata=treat$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 5.756853 0.3024978 6.649871 0.001
Residual 46 13.274204 0.6975022 NA NA
Total 49 19.031057 1.0000000 NA NA
6.3.4.2.4 Phylogenetic
betadisper(beta_div_phylo_treat$S, treat$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.58372 0.58372 35.413    999  0.001 ***
Residuals 48 0.79119 0.01648                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.001
2_Antibiotics    2.9795e-07              
adonis2(formula=beta_div_phylo_treat$S ~ time_point*Population, data=treat[labels(beta_div_phylo_treat$S),], permutations=999,strata=treat$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 2.947011 0.344846 8.070832 0.001
Residual 46 5.598866 0.655154 NA NA
Total 49 8.545877 1.000000 NA NA
6.3.4.2.5 Functional
betadisper(beta_div_func_treat$S, treat$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.17618 0.17618 4.7941    999  0.033 *
Residuals 48 1.76400 0.03675                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.026
2_Antibiotics      0.033451              
adonis2(formula=beta_div_func_treat$S ~ time_point*Population, data=treat[labels(beta_div_func_treat$S),], permutations=999,strata=treat$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 1.795938 0.3810423 9.439497 0.001
Residual 46 2.917286 0.6189577 NA NA
Total 49 4.713224 1.0000000 NA NA
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Does the FMT work?

6.3.5.1 Comparison between FMT2 vs Post-FMT2

#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant3<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
  column_to_rownames("newID")

transplant3_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

full_counts<-temp_genome_counts %>%
    t()%>%
    as.data.frame()%>%
    rownames_to_column("Tube_code")%>%
    full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant3_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))
6.3.5.1.1 Number of samples used
[1] 49
beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)
#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]
6.3.5.1.2 Richness
adonis2(formula=beta_div_richness_transplant3$S ~ Population+time_point+type, data=transplant3[labels(beta_div_richness_transplant3$S),], permutations=999, strata=transplant3$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.500812 0.2535872 5.096117 0.001
Residual 45 10.304350 0.7464128 NA NA
Total 48 13.805162 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_transplant3$S,transplant3_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 1.4169018 5.739828 0.15622903 0.001 0.003 *
Control vs Hot_control 1 2.0940966 8.509112 0.21005427 0.001 0.003 *
Treatment vs Hot_control 1 0.3004618 1.265034 0.04179854 0.177 0.531
6.3.5.1.3 Neutral
adonis2(formula=beta_div_neutral_transplant3$S ~ Population+time_point+type, data=transplant3[labels(beta_div_neutral_transplant3$S),], permutations=999, strata=transplant3$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 4.128749 0.3031142 6.524331 0.001
Residual 45 9.492350 0.6968858 NA NA
Total 48 13.621099 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_transplant3$S,transplant3_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 1.8758788 8.282671 0.21084796 0.001 0.003 *
Control vs Hot_control 1 2.4396317 10.635546 0.24945256 0.001 0.003 *
Treatment vs Hot_control 1 0.3158428 1.394345 0.04587515 0.122 0.366
6.3.5.1.4 Phylogenetic
adonis2(formula=beta_div_phylo_transplant3$S ~ Population+time_point+type, data=transplant3[labels(beta_div_phylo_transplant3$S),], permutations=999, strata=transplant3$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.3971179 0.2701357 5.551766 0.001
Residual 45 1.0729504 0.7298643 NA NA
Total 48 1.4700683 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_transplant3$S,transplant3_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.14387705 5.735321 0.15612552 0.001 0.003 *
Control vs Hot_control 1 0.22715701 9.044894 0.22036587 0.001 0.003 *
Treatment vs Hot_control 1 0.04648319 1.704277 0.05550617 0.125 0.375
6.3.5.1.5 Functional
adonis2(formula=beta_div_func_transplant3$S ~ Population+time_point+type, data=transplant3[labels(beta_div_func_transplant3$S),], permutations=999, strata=transplant3$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.0880056 0.0736928 1.193332 0.465
Residual 45 1.1062168 0.9263072 NA NA
Total 48 1.1942224 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_transplant3$S,transplant3_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.08177408 4.84137651 0.135077862 0.067 0.201
Control vs Hot_control 1 0.05266301 2.16167342 0.063277738 0.199 0.597
Treatment vs Hot_control 1 -0.00189892 -0.06088838 -0.002104017 0.845 1.000
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))
p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

6.3.5.2 Comparison between the different experimental time points (Acclimation vs Transplant samples)

The estimated time for calculating the 5151 pairwise combinations is 52 seconds.
ggarrange( p1, p2, p3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

6.3.5.3 Comparison of acclimation samples to transplant samples

transplant7<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1")%>%
  column_to_rownames("newID")

transplant7_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1")

transplant7_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

transplant7_counts <- transplant7_counts[, !names(transplant7_counts) %in% c("AD45 _ LI1_2nd_2", "AD48 _ LI1_2nd_6")]

identical(sort(colnames(transplant7_counts)),sort(as.character(rownames(transplant7))))
[1] TRUE
6.3.5.3.1 Number of samples used
[1] 73
beta_div_richness_transplant7<-hillpair(data=transplant7_counts, q=0)
beta_div_neutral_transplant7<-hillpair(data=transplant7_counts, q=1)
beta_div_phylo_transplant7<-hillpair(data=transplant7_counts, q=1, tree=genome_tree)
beta_div_func_transplant7<-hillpair(data=transplant7_counts, q=1, dist=dist)
#Arrange of metadata dataframe
transplant7_arrange<-transplant7[labels(beta_div_neutral_transplant7$S),]
transplant7_arrange <- transplant7_arrange %>%
  mutate(time_point = recode(time_point,
                             "3_Transplant1" = "Transplant",
                             "4_Transplant2" = "Transplant"))

transplant7_arrange$type_time <- interaction(transplant7_arrange$type, transplant7_arrange$time_point)
6.3.5.3.2 Richness
adonis2(formula=beta_div_richness_transplant7$S ~ Population*time_point+type, data=transplant7[labels(beta_div_richness_transplant7$S),], permutations=999, strata=transplant7$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 6 5.309519 0.2518733 3.703392 0.001
Residual 66 15.770599 0.7481267 NA NA
Total 72 21.080119 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_transplant7$S,transplant7_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.36208146 1.0521088 0.06169963 0.306 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.28008774 4.6054436 0.22350616 0.001 0.015 .
Control.1_Acclimation vs Control.Transplant 1 0.55038651 2.2107376 0.08124505 0.007 0.105
Control.1_Acclimation vs Treatment.Transplant 1 1.62289430 6.7106689 0.25123553 0.001 0.015 .
Control.1_Acclimation vs Hot_control.Transplant 1 1.73215888 7.4315069 0.25250175 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.36066298 5.0871520 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Control.Transplant 1 0.52860586 2.1820402 0.08027507 0.002 0.030 .
Treatment.1_Acclimation vs Treatment.Transplant 1 1.76810026 7.5736721 0.27467042 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.Transplant 1 1.87790626 8.3291875 0.27462613 0.001 0.015 .
Hot_control.1_Acclimation vs Control.Transplant 1 1.75314247 8.7706781 0.25971282 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.Transplant 1 0.27700454 1.5346880 0.07126586 0.067 1.000
Hot_control.1_Acclimation vs Hot_control.Transplant 1 0.26448976 1.4916174 0.06349573 0.095 1.000
Control.Transplant vs Treatment.Transplant 1 2.30884687 12.4299510 0.30002331 0.001 0.015 .
Control.Transplant vs Hot_control.Transplant 1 2.50396161 13.6713271 0.30604256 0.001 0.015 .
Treatment.Transplant vs Hot_control.Transplant 1 0.01688622 0.1023282 0.00392027 1.000 1.000
6.3.5.3.3 Neutral
adonis2(formula=beta_div_neutral_transplant7$S ~ Population+time_point*type, data=transplant7[labels(beta_div_neutral_transplant7$S),], permutations=999, strata=transplant7$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 8 7.284378 0.3492417 4.293351 0.001
Residual 64 13.573319 0.6507583 NA NA
Total 72 20.857698 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_transplant7$S,transplant7_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.23160196 0.7712905 0.045988741 0.694 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.40153474 5.7562378 0.264578733 0.001 0.015 .
Control.1_Acclimation vs Control.Transplant 1 0.56111203 2.5583085 0.092832565 0.002 0.030 .
Control.1_Acclimation vs Treatment.Transplant 1 1.88709838 8.3257794 0.293929402 0.001 0.015 .
Control.1_Acclimation vs Hot_control.Transplant 1 2.02585000 9.2317432 0.295588471 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.63477039 6.8326887 0.299250291 0.001 0.015 .
Treatment.1_Acclimation vs Control.Transplant 1 0.61335323 2.8313912 0.101733730 0.002 0.030 .
Treatment.1_Acclimation vs Treatment.Transplant 1 2.10939140 9.4473664 0.320822116 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.Transplant 1 2.24827218 10.3907678 0.320794118 0.001 0.015 .
Hot_control.1_Acclimation vs Control.Transplant 1 1.87351542 10.3925002 0.293635661 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.Transplant 1 0.34276062 1.9273510 0.087897118 0.044 0.660
Hot_control.1_Acclimation vs Hot_control.Transplant 1 0.31638309 1.8072337 0.075911118 0.055 0.825
Control.Transplant vs Treatment.Transplant 1 2.48701901 14.0199769 0.325894571 0.001 0.015 .
Control.Transplant vs Hot_control.Transplant 1 2.75304261 15.6912860 0.336064549 0.001 0.015 .
Treatment.Transplant vs Hot_control.Transplant 1 0.01764676 0.1022118 0.003915827 0.996 1.000
6.3.5.3.4 Phylogenetic
adonis2(formula=beta_div_phylo_transplant7$S ~ Population+time_point+type, data=transplant7[labels(beta_div_phylo_transplant7$S),], permutations=999, strata=transplant7$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 4 0.7377029 0.1879202 3.933904 0.022
Residual 68 3.1879143 0.8120798 NA NA
Total 72 3.9256172 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_transplant7$S,transplant7_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.04186923 0.43916424 0.026714511 0.722 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.15609416 2.55468892 0.137684276 0.030 0.450
Control.1_Acclimation vs Control.Transplant 1 0.03888650 0.83961027 0.032493148 0.487 1.000
Control.1_Acclimation vs Treatment.Transplant 1 0.28946588 4.58406811 0.186464994 0.002 0.030 .
Control.1_Acclimation vs Hot_control.Transplant 1 0.31864880 5.37781508 0.196429666 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.05218385 0.202081922 0.002 0.030 .
Treatment.1_Acclimation vs Control.Transplant 1 0.11794420 2.69844074 0.097422117 0.045 0.675
Treatment.1_Acclimation vs Treatment.Transplant 1 0.37640156 6.28511923 0.239113210 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.Transplant 1 0.40433696 7.18306079 0.246138020 0.001 0.015 .
Hot_control.1_Acclimation vs Control.Transplant 1 0.11597038 5.32063275 0.175478948 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.Transplant 1 0.03673004 1.13023077 0.053488804 0.363 1.000
Hot_control.1_Acclimation vs Hot_control.Transplant 1 0.04097680 1.30539166 0.056012432 0.264 1.000
Control.Transplant vs Treatment.Transplant 1 0.21736741 7.59281199 0.207494630 0.001 0.015 .
Control.Transplant vs Hot_control.Transplant 1 0.25837791 9.19762187 0.228810100 0.001 0.015 .
Treatment.Transplant vs Hot_control.Transplant 1 0.00180330 0.04804393 0.001844435 0.969 1.000
6.3.5.3.5 Functional
adonis2(formula=beta_div_func_transplant7$S ~ Population+time_point+type, data=transplant7[labels(beta_div_func_transplant7$S),], permutations=999, strata=transplant7$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 4 0.5122014 0.2427344 5.449191 0.038
Residual 68 1.5979298 0.7572656 NA NA
Total 72 2.1101312 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_transplant7$S,transplant7_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.0905830704 1.66462866 0.0942351347 0.247 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.0861813922 1.63467278 0.0926965190 0.238 1.000
Control.1_Acclimation vs Control.Transplant 1 0.0706284677 2.02459114 0.0749166241 0.178 1.000
Control.1_Acclimation vs Treatment.Transplant 1 0.3227173802 7.20965350 0.2649667516 0.005 0.075
Control.1_Acclimation vs Hot_control.Transplant 1 0.3449345536 8.46661273 0.2778980651 0.005 0.075
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.0010225901 0.05430389 0.0033825127 0.658 1.000
Treatment.1_Acclimation vs Control.Transplant 1 -0.0046542303 -0.35270812 -0.0143102181 0.801 1.000
Treatment.1_Acclimation vs Treatment.Transplant 1 0.0783171063 4.43726923 0.1815779491 0.073 1.000
Treatment.1_Acclimation vs Hot_control.Transplant 1 0.0836921311 5.20043693 0.1911894629 0.060 0.900
Hot_control.1_Acclimation vs Control.Transplant 1 -0.0042700245 -0.35258632 -0.0143052054 0.813 1.000
Hot_control.1_Acclimation vs Treatment.Transplant 1 0.0824858621 5.06251874 0.2019956092 0.063 0.945
Hot_control.1_Acclimation vs Hot_control.Transplant 1 0.0887346857 5.97129920 0.2134795084 0.049 0.735
Control.Transplant vs Treatment.Transplant 1 0.1927489878 15.76935832 0.3522355226 0.004 0.060
Control.Transplant vs Hot_control.Transplant 1 0.2075592800 18.09824701 0.3686128958 0.003 0.045 .
Treatment.Transplant vs Hot_control.Transplant 1 -0.0001900114 -0.01304792 -0.0005020952 0.702 1.000
beta_richness_nmds_transplant7 <- beta_div_richness_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant7 <- beta_div_neutral_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant7 <- beta_div_phylo_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant7 <- beta_div_func_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))
p0<-beta_richness_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p1<-beta_neutral_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p2<-beta_phylo_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p3<-beta_func_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
  theme_classic()+
  theme(legend.position="none")

6.3.5.4 Comparison between Acclimation vs Post-FMT1

post3 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post3))]
identical(sort(colnames(post3.counts)),sort(as.character(rownames(post3))))

post3_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")
6.3.5.4.1 Number of samples used
[1] 53
beta_div_richness_post3<-hillpair(data=post3.counts, q=0)
beta_div_neutral_post3<-hillpair(data=post3.counts, q=1)
beta_div_phylo_post3<-hillpair(data=post3.counts, q=1, tree=genome_tree)
beta_div_func_post3<-hillpair(data=post3.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post3_arrange<-post3[labels(beta_div_neutral_post3$S),]
post3_arrange$type_time <- interaction(post3_arrange$type, post3_arrange$time_point)
6.3.5.4.2 Richness
betadisper(beta_div_richness_post3$S, post3_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     2 0.099607 0.049803 9.5441    999  0.001 ***
Residuals 50 0.260911 0.005218                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Control Hot_control Treatment
Control                 0.00100000     0.894
Hot_control 0.00102653                 0.001
Treatment   0.88832670  0.00010131          
adonis2(formula=beta_div_richness_post3$S ~ time_point*type, data=post3[labels(beta_div_neutral_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 4.403362 0.2369995 2.919782 0.001
Residual 47 14.176268 0.7630005 NA NA
Total 52 18.579631 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_richness_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3620815 1.052109 0.06169963 0.339 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2800877 4.605444 0.22350616 0.001 0.015 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.6845657 1.998114 0.11101796 0.005 0.075
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8437461 2.499232 0.14281954 0.001 0.015 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.1208022 3.568670 0.18236649 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.7216200 2.172734 0.11956009 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9551308 2.926054 0.16322910 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.2263345 4.039487 0.20157637 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.4319792 5.384836 0.25180628 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8172413 3.194690 0.17558364 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.5796135 2.441615 0.13239702 0.002 0.030 .
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.016 0.240
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.121 1.000
6.3.5.4.3 Neutral
betadisper(beta_div_neutral_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00945 0.0094472 1.1428    999  0.267
Residuals 51 0.42161 0.0082669                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                      0.26
5_Post-FMT1          0.2901            
adonis2(formula=beta_div_neutral_post3$S ~ time_point*type, data=post3[labels(beta_div_neutral_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 5.302354 0.3027004 4.080576 0.001
Residual 47 12.214484 0.6972996 NA NA
Total 52 17.516838 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_neutral_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2316020 0.7712905 0.04598874 0.715 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4015347 5.7562378 0.26457873 0.001 0.015 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.8332162 2.9081103 0.15380227 0.001 0.015 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 1.1719595 4.0685514 0.21336447 0.001 0.015 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.4260875 5.2413171 0.24675104 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.9517634 3.3715700 0.17404733 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 1.3127773 4.6298256 0.23585668 0.002 0.030 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.6713369 6.2395460 0.28056085 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.5409781 6.8338056 0.29928456 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9133614 4.0964534 0.21451383 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.6954835 3.2951234 0.17077493 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.2508491 0.13047758 0.019 0.285
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.1436369 0.20570451 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.6372683 0.09840968 0.050 0.750
6.3.5.4.4 Phylogenetic
betadisper(beta_div_phylo_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.05132 0.051320 2.6745    999  0.108
Residuals 51 0.97861 0.019189                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.099
5_Post-FMT1         0.10812            
adonis2(formula=beta_div_phylo_post3$S ~ time_point*type, data=post3[labels(beta_div_phylo_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 0.7935087 0.2278753 2.774199 0.006
Residual 47 2.6886978 0.7721247 NA NA
Total 52 3.4822065 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_phylo_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.04186923 0.4391642 0.02671451 0.749 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.15609416 2.5546889 0.13768428 0.033 0.495
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.19193367 2.9749922 0.15678490 0.019 0.285
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.14627288 1.7907381 0.10665035 0.164 1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.25061348 3.6146185 0.18428187 0.011 0.165
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.0521838 0.20208192 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.26358465 4.3608960 0.21417997 0.003 0.045 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.25319427 3.2738422 0.17915456 0.043 0.645
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.39050120 5.9837393 0.27218933 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 0.14203376 5.4200212 0.25303529 0.002 0.030 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09666753 2.3682173 0.13635351 0.012 0.180
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.09252600 2.9824958 0.15711821 0.009 0.135
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.784 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.105 1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.733 1.000
6.3.5.4.5 Functional
betadisper(beta_div_func_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00554 0.0055401 0.2063    999  0.629
Residuals 51 1.36938 0.0268505                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.628
5_Post-FMT1         0.65159            
adonis2(formula=beta_div_func_post3$S ~ time_point*type, data=post3[labels(beta_div_func_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 0.2935075 0.1510442 1.672425 0.054
Residual 47 1.6496826 0.8489558 NA NA
Total 52 1.9431901 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_func_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.090583070 1.66462866 0.094235135 0.197 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.086181392 1.63467278 0.092696519 0.213 1.000
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.028641941 0.50417680 0.030548437 0.553 1.000
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.234795406 4.03037749 0.211786524 0.047 0.705
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.134726259 2.20299547 0.121023788 0.171 1.000
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.001022590 0.05430389 0.003382513 0.658 1.000
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.002157067 0.09411569 0.005847832 0.602 1.000
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.056602363 2.56037069 0.145803909 0.185 1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.009569124 0.35095521 0.021463896 0.502 1.000
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 -0.001745663 -0.08225018 -0.005167199 0.724 1.000
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.057758674 2.84545622 0.159449901 0.162 1.000
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.005575266 0.21803560 0.013444020 0.541 1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.119540855 4.84764704 0.244242909 0.075 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.052587837 1.77308932 0.099762584 0.246 1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.012980354 0.44307662 0.028690955 0.469 1.000
beta_richness_nmds_post3 <- beta_div_richness_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- beta_div_neutral_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- beta_div_phylo_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- beta_div_func_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

6.3.5.5 Comparison between Acclimation vs Post-FMT2

post4 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post4.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post4))]
identical(sort(colnames(post4.counts)),sort(as.character(rownames(post4))))

post4_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")
6.3.5.5.1 Number of samples used
[1] 54
beta_div_richness_post4<-hillpair(data=post4.counts, q=0)
beta_div_neutral_post4<-hillpair(data=post4.counts, q=1)
beta_div_phylo_post4<-hillpair(data=post4.counts, q=1, tree=genome_tree)
beta_div_func_post4<-hillpair(data=post4.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post4_arrange<-post4[labels(beta_div_neutral_post4$S),]
post4_arrange$type_time <- interaction(post4_arrange$type, post4_arrange$time_point)
6.3.5.5.2 Richness
betadisper(beta_div_richness_post4$S, post4_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.06809 0.034047 3.8471    999  0.031 *
Residuals 51 0.45135 0.008850                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0360000     0.891
Hot_control 0.0349385                 0.003
Treatment   0.8855174   0.0047257          
adonis2(formula=beta_div_richness_post4$S ~ time_point*Population, data=post4[labels(beta_div_neutral_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.310172 0.1883377 3.867324 0.001
Residual 50 14.265560 0.8116623 NA NA
Total 53 17.575732 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_richness_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3620815 1.052109 0.06169963 0.305 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2800877 4.605444 0.22350616 0.001 0.015 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.8430295 2.845779 0.15100353 0.001 0.015 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.5232174 1.683240 0.09518843 0.022 0.330
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.1217138 3.634271 0.18509835 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.9130048 3.195028 0.16645080 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.5959230 1.984036 0.11032208 0.002 0.030 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.2747787 4.275366 0.21086503 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6397330 2.913695 0.15405213 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.4575447 6.224524 0.28007456 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3276169 1.412318 0.08111028 0.040 0.600
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.002 0.030 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.015 .
6.3.5.5.3 Neutral
betadisper(beta_div_neutral_post4$S, post4_arrange$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01544 0.0154447 2.0972    999  0.163
Residuals 52 0.38294 0.0073643                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.161
6_Post-FMT2         0.15357            
adonis2(formula=beta_div_neutral_post4$S ~ time_point*Population, data=post4[labels(beta_div_neutral_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.863228 0.229321 4.959284 0.001
Residual 50 12.983151 0.770679 NA NA
Total 53 16.846379 1.000000 NA NA
pairwise <- pairwise.adonis(beta_div_neutral_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2316020 0.7712905 0.04598874 0.724 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4015347 5.7562378 0.26457873 0.001 0.015 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 1.1746426 4.5564741 0.22165640 0.001 0.015 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.5286441 1.9819408 0.11021840 0.005 0.075
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.3443224 4.9104417 0.23483204 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 1.3540292 5.3398081 0.25022756 0.002 0.030 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.6311089 2.4041625 0.13063146 0.002 0.030 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.6125755 5.9825981 0.27215155 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6202327 3.1519868 0.16457754 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.5701179 7.6327037 0.32297209 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3634438 1.7083388 0.09647087 0.031 0.465
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.6483346 0.22511910 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.2065321 0.12119453 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.7710313 0.26507845 0.001 0.015 .
6.3.5.5.4 Phylogenetic
betadisper(beta_div_phylo_post4$S, post4$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.06978 0.069777 5.0345    999  0.028 *
Residuals 52 0.72071 0.013860                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.024
6_Post-FMT2        0.029131            
adonis2(formula=beta_div_phylo_post4$S ~ time_point*Population, data=post4[labels(beta_div_phylo_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.757493 0.2376349 5.195124 0.001
Residual 50 2.430141 0.7623651 NA NA
Total 53 3.187634 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_phylo_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.04186923 0.4391642 0.02671451 0.735 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.15609416 2.5546889 0.13768428 0.038 0.570
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.26322331 4.3060281 0.21205664 0.007 0.105
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.16047895 2.5405742 0.13702781 0.037 0.555
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.25529510 4.0109138 0.20043631 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.0521838 0.20208192 0.005 0.075
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.36496892 6.3966666 0.28560797 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.22628210 3.8292220 0.19311005 0.014 0.210
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.34830814 5.8463335 0.26761166 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.10002871 4.3836237 0.21505615 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 0.12577510 5.0601287 0.24027055 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.06334378 2.4997737 0.13512455 0.018 0.270
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.3820253 0.12958449 0.027 0.405
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.7224602 0.14541146 0.007 0.105
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.0436561 0.20174244 0.002 0.030 .
6.3.5.5.5 Functional
betadisper(beta_div_func_post4$S, post4$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00527 0.005269 0.1889    999  0.665
Residuals 52 1.45058 0.027896                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.657
6_Post-FMT2         0.66565            
adonis2(formula=beta_div_func_post4$S ~ time_point*Population, data=post4[labels(beta_div_func_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.0773959 0.0417692 0.726498 0.31
Residual 50 1.7755477 0.9582308 NA NA
Total 53 1.8529436 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_func_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.0905830704 1.664628661 0.0942351347 0.222 1
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.0861813922 1.634672780 0.0926965190 0.253 1
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.1197900330 2.213130846 0.1215129274 0.164 1
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.1125623700 2.150784454 0.1184953995 0.179 1
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0657004998 0.954588109 0.0563026423 0.284 1
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.0010225901 0.054303886 0.0033825127 0.634 1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 -0.0005177706 -0.025585400 -0.0016016487 0.710 1
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.0013301207 0.072110871 0.0044867082 0.609 1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0060959077 0.174487757 0.0107878382 0.575 1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.0010345754 0.055797964 0.0034752533 0.654 1
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 -0.0001056284 -0.006306177 -0.0003942915 0.679 1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0017235602 0.051851181 0.0032302306 0.744 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -0.0080428882 -0.442986255 -0.0284750185 0.847 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -0.0011796256 -0.034047378 -0.0021324990 0.902 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.0036300838 0.110487573 0.0068581148 0.673 1
beta_richness_nmds_post4 <- beta_div_richness_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- beta_div_neutral_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- beta_div_phylo_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- beta_div_func_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

6.3.6 6. Are there differences between the control and the treatment group?

6.3.6.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")
6.3.6.1.1 Number of samples used
[1] 26
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]
6.3.6.1.2 Richness
betadisper(beta_div_richness_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.017675 0.0088373 2.3825    999  0.111
Residuals 23 0.085312 0.0037092                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0050000     0.682
Hot_control 0.0068795                 0.229
Treatment   0.6248469   0.2084296          
adonis2(formula=beta_div_richness_post1$S ~ Population+type, data=post1[labels(beta_div_richness_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.195567 0.1448246 1.947534 0.001
Residual 23 7.059710 0.8551754 NA NA
Total 25 8.255277 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.013      0.039   .
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.106      0.318    
6.3.6.1.3 Neutral
betadisper(beta_div_neutral_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.575
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.23900     0.966
Hot_control 0.21166                 0.451
Treatment   0.95468     0.43604          
adonis2(formula=beta_div_neutral_post1$S ~ Population+type, data=post1[labels(beta_div_neutral_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.395968 0.1900228 2.697931 0.001
Residual 23 5.950350 0.8099772 NA NA
Total 25 7.346318 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.010      0.030   .
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.002      0.006   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.059      0.177    
6.3.6.1.4 Phylogenetic
betadisper(beta_div_phylo_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999    0.9
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92000     0.702
Hot_control 0.91505                 0.758
Treatment   0.63312     0.73046          
adonis2(formula=beta_div_phylo_post1$S ~ Population+type, data=post1[labels(beta_div_phylo_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.0745104 0.0705947 0.8735033 0.532
Residual 23 0.9809570 0.9294053 NA NA
Total 25 1.0554673 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.798      1.000    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.119      0.357    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.681      1.000    
6.3.6.1.5 Functional
betadisper(beta_div_func_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00400 0.0019999 0.1431    999  0.862
Residuals 23 0.32135 0.0139717                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.61000     0.707
Hot_control 0.60188                 0.848
Treatment   0.74597     0.84473          
adonis2(formula=beta_div_func_post1$S ~ Population+type, data=post1[labels(beta_div_func_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.1230554 0.1608583 2.204479 0.16
Residual 23 0.6419374 0.8391417 NA NA
Total 25 0.7649929 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11954085 4.8476470 0.24424291   0.070      0.210    
2   Control vs Hot_control  1 0.05258784 1.7730893 0.09976258   0.232      0.696    
3 Treatment vs Hot_control  1 0.01298035 0.4430766 0.02869096   0.473      1.000    
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.6.2 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")
6.3.6.2.1 Number of samples used
[1] 27
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]
6.3.6.2.2 Richness
betadisper(beta_div_richness_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.832
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.69800     0.814
Hot_control 0.67789                 0.649
Treatment   0.79246     0.59820          
adonis2(formula=beta_div_richness_post2$S ~ type, data=post2[labels(beta_div_richness_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.504341 0.1967776 2.939822 0.001
Residual 24 6.140538 0.8032224 NA NA
Total 26 7.644879 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 0.6463814 2.560441 0.1379515 0.001 0.003 *
Treatment vs Hot_control 1 0.4796256 1.916520 0.1069694 0.001 0.003 *
Control vs Hot_control 1 1.1305044 4.268317 0.2105906 0.001 0.003 *
6.3.6.2.3 Neutral
betadisper(beta_div_neutral_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.492
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.47700     0.655
Hot_control 0.44675                 0.267
Treatment   0.65989     0.25095          
adonis2(formula=beta_div_neutral_post2$S ~ type, data=post2[labels(beta_div_neutral_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.923807 0.2603795 4.224537 0.001
Residual 24 5.464666 0.7396205 NA NA
Total 26 7.388473 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 1.0227481 4.648335 0.2251191 0.001 0.003 *
Treatment vs Hot_control 1 0.5010202 2.206532 0.1211945 0.002 0.006 *
Control vs Hot_control 1 1.3619424 5.771031 0.2650785 0.001 0.003 *
6.3.6.2.4 Phylogenetic
betadisper(beta_div_phylo_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.954
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93000     0.849
Hot_control 0.93765                 0.761
Treatment   0.83933     0.76015          
adonis2(formula=beta_div_phylo_post2$S ~ type, data=post2[labels(beta_div_phylo_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.1594363 0.2042241 3.079623 0.001
Residual 24 0.6212564 0.7957759 NA NA
Total 26 0.7806927 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 0.05927454 2.382025 0.1295845 0.034 0.102
Treatment vs Hot_control 1 0.06906280 2.722460 0.1454115 0.004 0.012 .
Control vs Hot_control 1 0.11081709 4.043656 0.2017424 0.001 0.003 *
6.3.6.2.5 Functional
betadisper(beta_div_func_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01259 0.0062962 0.3249    999  0.798
Residuals 24 0.46507 0.0193778                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.54400     0.653
Hot_control 0.45381                 0.798
Treatment   0.57452     0.74365          
adonis2(formula=beta_div_func_post2$S ~ type, data=post2[labels(beta_div_func_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 -0.0037283 -0.0054704 -0.065288 0.915
Residual 24 0.6852623 1.0054704 NA NA
Total 26 0.6815340 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 -0.008042888 -0.44298625 -0.028475019 0.842 1
Treatment vs Hot_control 1 -0.001179626 -0.03404738 -0.002132499 0.902 1
Control vs Hot_control 1 0.003630084 0.11048757 0.006858115 0.681 1
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.6.3 Post1 vs Post2

post5 <- meta %>%
  filter(time_point == "6_Post-FMT2" | time_point == "5_Post-FMT1")

post5.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post5))]
identical(sort(colnames(post5.counts)),sort(as.character(rownames(post5))))

post5_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")
6.3.6.3.1 Number of samples used
[1] 53
beta_div_richness_post5<-hillpair(data=post5.counts, q=0)
beta_div_neutral_post5<-hillpair(data=post5.counts, q=1)
beta_div_phylo_post5<-hillpair(data=post5.counts, q=1, tree=genome_tree)
beta_div_func_post5<-hillpair(data=post5.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post5_arrange<-post5[labels(beta_div_neutral_post5$S),]
post5_arrange$type_time <- interaction(post5_arrange$type, post5_arrange$time_point)
6.3.6.3.2 Richness
betadisper(beta_div_richness_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01841 0.0092048 1.7364    999  0.187
Residuals 50 0.26505 0.0053010                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.045000     0.699
Hot_control 0.039117                 0.202
Treatment   0.716358    0.218648          
adonis2(formula=beta_div_richness_post5$S ~ type, data=post5[labels(beta_div_richness_post5$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 2.28826 0.1390012 4.036044 0.001
Residual 50 14.17390 0.8609988 NA NA
Total 52 16.46216 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.020 0.300
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.7628135 2.683925 0.14364890 0.002 0.030 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3432605 1.148733 0.06698647 0.258 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.1269580 3.799256 0.19188884 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.118 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3571397 1.297184 0.07959561 0.119 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.7769467 2.670898 0.15114670 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.6502360 2.253407 0.13060650 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4132091 1.616138 0.09174188 0.009 0.135
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0163992 3.760571 0.19030682 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.2732563 1.019281 0.05988979 0.436 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.015 .
6.3.6.3.3 Neutral
betadisper(beta_div_neutral_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.01992 0.0099587 1.565    999  0.233
Residuals 50 0.31818 0.0063636                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.12300     0.868
Hot_control 0.10701                 0.162
Treatment   0.87156     0.17449          
adonis2(formula=beta_div_neutral_post5$S ~ type, data=post5[labels(beta_div_neutral_post5$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 2.928527 0.1906221 5.887921 0.001
Residual 50 12.434468 0.8093779 NA NA
Total 52 15.362995 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.250849 0.13047758 0.015 0.225
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.143637 0.20570451 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.8908158 3.714692 0.18842252 0.001 0.015 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3860927 1.552176 0.08843210 0.082 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.3122237 5.130273 0.24279254 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.637268 0.09840968 0.039 0.585
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3157079 1.325203 0.08117526 0.175 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0579520 4.270010 0.22158835 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.7454015 2.920049 0.16294873 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4377161 1.942126 0.10824392 0.009 0.135
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.3766597 5.875279 0.26858075 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.3176516 1.316137 0.07600637 0.167 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.648335 0.22511910 0.002 0.030 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.206532 0.12119453 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.771031 0.26507845 0.001 0.015 .
6.3.6.3.4 Phylogenetic
betadisper(beta_div_phylo_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00051 0.0002543 0.0265    999  0.972
Residuals 50 0.47996 0.0095993                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.89200     0.831
Hot_control 0.88926                 0.926
Treatment   0.82391     0.91902          
adonis2(formula=beta_div_phylo_post5$S ~ type, data=post5[labels(beta_div_phylo_post5$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.1778594 0.0910216 2.503403 0.007
Residual 50 1.7761762 0.9089784 NA NA
Total 52 1.9540356 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.779 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.111 1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07917244 3.0180046 0.15869197 0.009 0.135
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.04335491 1.5335604 0.08746429 0.190 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.10783045 3.7500438 0.18987521 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.681 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.06393539 1.5651817 0.09448624 0.147 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.05265949 1.2240203 0.07544494 0.291 1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.09753501 2.2402429 0.12994265 0.018 0.270
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07228545 2.3279593 0.12701683 0.036 0.540
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.11759094 3.5538444 0.18174658 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.06667255 1.9859527 0.11041687 0.100 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.3820253 0.12958449 0.020 0.300
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.7224602 0.14541146 0.005 0.075
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.0436561 0.20174244 0.001 0.015 .
6.3.6.3.5 Functional
betadisper(beta_div_func_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00785 0.0039232 0.2322    999  0.791
Residuals 50 0.84483 0.0168966                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.53800     0.605
Hot_control 0.52384                 0.849
Treatment   0.58787     0.85068          
adonis2(formula=beta_div_func_post5$S ~ type, data=post5[labels(beta_div_func_post5$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.0668661 0.0466009 1.221967 0.315
Residual 50 1.3680018 0.9533991 NA NA
Total 52 1.4348679 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.1195408549 4.84764704 0.2442429086 0.072 1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.0525878365 1.77308932 0.0997625840 0.245 1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.0265995825 1.17541806 0.0684360667 0.293 1
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0145818992 0.69975992 0.0419023938 0.417 1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -0.0080695208 -0.21617323 -0.0136958691 0.914 1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.0129803540 0.44307662 0.0286909552 0.467 1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.0267162134 1.22560581 0.0755352882 0.321 1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0384388433 1.93281582 0.1141461550 0.217 1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.0553988290 1.47819391 0.0897060633 0.235 1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 -0.0040061386 -0.14850469 -0.0093684974 0.767 1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0024023972 0.09538980 0.0059265296 0.627 1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -0.0004960759 -0.01190328 -0.0007445087 0.838 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -0.0080428882 -0.44298625 -0.0284750185 0.846 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -0.0011796256 -0.03404738 -0.0021324990 0.921 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.0036300838 0.11048757 0.0068581148 0.697 1
beta_richness_nmds_post5 <- beta_div_richness_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post5 <- beta_div_neutral_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post5 <- beta_div_phylo_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post5 <- beta_div_func_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")